% estimateIRF 
% this version: march 2010 (draws on mkimp_restricted)

function [impzout,erzout,azeroout,a0betazout,var_explode]...
    = estimateIRF(vardata,nlags,nstep,teststationarity,exo);

[nobs,nvars] = size(vardata);  

y   = vardata((nlags+1):end,1:end);
x   = zeros(length(y),nlags*nvars);
erz = zeros(size(y));
exo = exo((nlags+1):end,:);

for zx = 1:nlags;
     x(:,(zx-1)*nvars+1:zx*nvars) = vardata((nlags+1)-zx:end-zx,:);
end;

azero = eye(nvars);
betaz = zeros(nvars,nvars*nlags);
nexo = size(exo,2);
bexo=zeros(nvars,nexo);

for zx=1:nvars;
    xmat = [exo y(:,1:(zx-1)) x];
    btmp = xmat\y(:,zx);
    erz(:,zx) = y(:,zx)-xmat*btmp;
    bexo(zx,:)=btmp(1:nexo)';      
    btmp = btmp(1+nexo:end)';
    azero(zx,1:(zx-1)) = -1*btmp(1:(zx-1));
    btmp = btmp((zx-1)+1:end);
    betaz(zx,:) = btmp;
end

a0inv = inv(azero);
zlag=zeros(1,nlags*nvars);

a0betaz = zeros(size(betaz));
for mx=1:nlags;
   a0betaz(:,1+(mx-1)*nvars:mx*nvars) = a0inv*betaz(:,1+(mx-1)*nvars:mx*nvars);
end;

betaz= a0betaz';
a0betazout = [a0inv*bexo  a0betaz];
azeroout = azero;
erzout   = erz;
 
bigA = [betaz'; kron(eye(nvars),eye(nlags-1)) zeros(nvars*(nlags-1),nvars)];
var_explode = sum(abs(eig(bigA))>1);
if teststationarity==1;
if var_explode > 0;
  disp('Warning: VAR Explosive')
elseif var_explode <=0
  disp('VAR stationary');
end
end


%calculate moving average representation and through this the impulse responses.
%in calculating moving average representation include the std of the shocks.

marep = zeros([nvars nvars nstep]);

for allshks = 1:size(erz,2); % # of shocks...
   %shock = [std(erz(:,allshks))*a0inv(:,allshks)';zeros(nstep-1,nvars)];    
   shock = [a0inv(:,allshks)';zeros(nstep-1,nvars)];  
   zlag=zeros(1,nlags*nvars);   
   for jj=1:nstep;
      imptmp = shock(jj,:)+zlag*betaz(1:length(betaz),:);
      marep(:,allshks,jj)=imptmp;
      zlag=[imptmp zlag(1,1:(nlags-1)*nvars)];
   end;
end;

for zx= 1:size(marep,2)
    impzout(:,:,zx) = squeeze(marep(:,zx,:))';
end;
